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Using a code that employs a self-consistent method for computing the effects of photoion¬ 
ization on circumstellar gas dynamics, we model the formation of wind-driven nebulae around 
massive Wolf-Rayet (W-R) stars. Our algorithm incorporates a simplified model of the photo¬ 
ionization source, computes the fractional ionization of hydrogen due to the photoionizing fiux 
and recombination, and determines self-consistently the energy balance due to ionization, photo¬ 
heating and radiative cooling. We take into account changes in stellar properties and mass-loss 
over the star’s evolution. Our multi-dimensional simulations clearly reveal the presence of 
strong ionization front instabilities. Using various X-ray emission models, and abundances con¬ 
sistent with those derived for W-R nebulae, we compute the X-ray fiux and spectra from our 
wind bubble models. We show the evolution of the X-ray spectral features with time over the 
evolution of the star, taking the absorption of the X-rays by the ionized bubble into account. 
Our simulated X-ray spectra compare reasonably well with observed spectra of Wolf-Rayet 
bubbles. They suggest that X-ray nebulae around massive stars may not be easily detectable, 
consistent with observations. 


1 Introduction: 

Massive stars (> 8 Mq) lose mass throughout their 
lifetime, via winds and eruptions, before ending their 
lives in a cataclysmic supernova (SN) explosion. The 
interaction of this material with the surrounding 
medium creates vast wind-blown cavities surrounded 
by a dense shell, referred to as wind-blown bub¬ 
bles. As the star evolves through various stages, 
the mass-loss parameters will change, affecting the 
structure of the bubble. When the star finally ex¬ 
plodes, the resulting SN shock wave will expand 
within the bubble, and the dynamics and kinemat¬ 
ics of the shock wave will depend on the bubble 
parameters (Dwarkadas (2005, 2007a,b)). Similarly, 
the relativistic blast waves associated with gamma- 
ray bursts are expected to expand within wind bub¬ 
bles surrounding Wolf-Rayet (W-R) stars. Using an 
ionization-gasdynamics code, AVATAR, we compute 
the structure and evolution of the wind-blown bub¬ 
bles around massive stars. Using the ISIS package, 
we compute the X-ray spectrum from the simula¬ 
tions as would be observed by the Chandra satellite, 
and compare to observational data. 


2 The AVATAR Code: 

We have further developed a code that combines 
photoionization from the star with the gasdynamics, 
and used it to compute the structure of wind bubbles 
around massive stars. The method operator splits 
the contribution due to photoionization effects from 
the usual gas dynamics, and utilizes a backward- 
Euler scheme together with a Newton-Raphson iter¬ 
ation procedure for achieving a solution. The effects 
of geometrical dilution and of column absorption of 


radiation are considered. The gasdynamic algorithm 
makes use of a multidimensional covariant imple¬ 
mentation of well established Eulerian finite differ¬ 
ence algorithms (Rosenberg (1995)). A second-order 
(van Leer) monotonic transport algorithm is used 
for advection of total mass and the neutral compo¬ 
nent, and a third order piecewise parabolic scheme 
is available. Tabulated functions are used to com¬ 
pute the collisional ionization rate and cooling func¬ 
tion. Shocks are treated using an artificial viscosity. 
Grid expansion is available to study fiow over dis¬ 
tance scales spanning several orders of magnitude. 
The algorithm incorporates a simplified model of the 
photo-ionization source, computes the fractional ion¬ 
ization of hydrogen due to the photo-ionizing fiux 
and recombination, and determines self-consistently 
the energy balance due to ionization, photo-heating 
and radiative cooling. In this, our method is superior 
to that of other calculations, such as Garcia-Segura 
et al. (1996); van Marie et al. (2005), who use the 
on-the-spot approximation and do not take the re¬ 
combination time into account. It is comparable to 
the work of Toala & Arthur (2011). 

3 Evolution of the Wind-Blown 
bubble around a 40 Mq star: 

A model for a 40 Mq star using the AVATAR code 
is shown in Eigure 1. Stellar parameters are adapted 
from van Marie et al. (2005). The evolution of the 
star can be divided into 3 main phases. An inhomo¬ 
geneous pressure and density distribution develops 
in the main-sequence (MS) stage, accompanied by 
vorticity deposition near the inner shock. The in¬ 
clusion of photo-ionization results in the formation 
of a dense, lower temperature (^ 10^ K) region of 
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ionized material outside the wind bubble during the 
MS phase. The nebula is fully ionized, and the ion¬ 
ization front is trapped in the dense shell, which is 
unstable to ionization instabilities. In the red su¬ 
pergiant (RSG) stage the surface temperature of the 
star decreases considerably. Consequently the ioniz¬ 
ing radiation drops considerably and recombination 
reduces the ionization fraction to ^ 40%, although 
this goes up again in the W-R stage. The high- 
density RSG wind is followed by a higher momentum 
W-R wind which breaks up the RSG shell, scattering 
the dense material into a turbulent W-R nebula. 
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Fig. 2: Maps showing the area from which X-ray emis¬ 
sion arises from the bubble in Figure 1. An increase of 
the X-ray flux in the W-R stage is seen. 



Fig. 1: Density snapshots from calculation of wind bub¬ 
ble evolution around 40 Mq star (600 x 400 zones), com¬ 
puted with AVATAR with expanding grid (Dwarkadas V 
Rosenberg (2013)). Time increases from left to right and 
top to bottom, and is listed at the top of each plot. 


We can use the simulations to illustrate the re¬ 
gions from which the X-ray flux arises. We assume 
that thermal bremsstrahlung is the dominant con¬ 
tributor to the X-ray flux, although line emission 
will also play a role. The flux is then proportional 
to rig where rie is the density and T the tem¬ 
perature. We show maps of this quantity at various 
times in Figure 2. The normalization is arbitrary, 
however they are all normalized to the same scale. 
Only zones having temperature T > 10^ K are as¬ 
sumed to emit X-rays. The densest clumps do not 
contribute because their temperature is lower than 
10^ K. It is clear that the highest X-ray flux arises in 
the W-R stage, and that the highest levels of X-ray 
emission emanate from very small regions. 


4 THE X-RAY SPECTRUM OF 
THE BUBBLE: 

In order to calculate the X-ray emission accurately, 
we use the ISIS package (Houck & Denicola (2000)). 
The data are read into ISIS, the spectrum calcu¬ 
lated for every zone, taking the absorption outside 
that zone into account, and added together. A point 
source is assumed. We use the VMEKAL model 
in XSPEC to model the spectrum, since it is the 
one most commonly used to fit the observed X-ray 
spectra. Herein we present (Eigure 3) X-ray spectra 
as would be detected by the Chandra ACIS instru¬ 
ment for an object at 1.5kpc distance, integrated 
over 50,000s. Eoreground absorption of 2. x 10^^ 
cm“^ is added to absorption due to the neutral ma¬ 
terial surrounding the bubble, which varies, and can 
be as high as 2. x 10^^ cm-2 (total NH in units of 
10^^ cm“^). Line broadening is based on the un¬ 
derlying fluid velocity. Solar abundances (Anders 
& Grevesse (1989)) are used for the MS and RSG 
stages. Abundances in the W-R phase are from Chu 
et al. (2003) for the W-R bubble S308. 

Our spectra are comparable in shape and temper¬ 
ature to observed spectra (Eigure 4). Our simula¬ 
tions show that the X-ray emission in the MS phase 
is generally too weak to be detected by current X- 
ray satellites. It is higher in the W-R stage, but 
still difficult to detect. The X-ray emission depends 
strongly on the density within the bubble, which is 
a function of the density of the material around the 
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Fig. 3: Simulated X-ray spectra for a point source calculated from the simulation shown in Figure 1. The leftmost 
image is from the MS phase, the next three are from the W-R phase. Note that as time increases in the W-R stage, 
the intensity of the spectrum (counts/sec/keV) also increases. 


Vmekal, 5308 Abund, NH = 0.055793466193511214 Time= 4.5828E+06 yns 
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Fig. 4: Simulated and Observed spectra. Panel 1 resembles observed Chandra spectrum for NGC 6888 (panel 3) 
and Panel 2 XMM spectrum of S308 (panel 4), from Chu et al. (2006). Although the spectral shapes are similar, 
the counts are significantly lower in our simulated spectra as compared to the data. This seems to be generally 
true over all our simulations, in agreement with the fact that diffuse X-ray emission is not generally observed. 


nebula. Thus it is possible that a higher ambient 
density leads to a higher intensity, although unless 
it was all ionized the absorption would correspond¬ 
ingly increase. The general low intensity is in agree¬ 
ment with observations, which have only detected 3 
W-R nebulae, and no MS nebulae around massive 
stars (Chu et al. (2006); Toala et al. (2014, 2015)) 
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